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The  paper  addresses  the  basic  notions  of  the  theory  of  feedback  and 
adaptive  finite  element  methods.  An  illustrative  example  of  the  finite 
element  method  for  solving  nonlinear  differential  equation  of  the  elliptic 
type  is  presented.  Numerical  example  of  regridding  process  is  given. 
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INTRODUCTION 


The  notion  of  feedback  and  adaptivity  in  numerical  computations  has 
attracted  great  attention  in  the  recent  years.  All  the  modern  ODE's  solvers 
for  Initial  value  problems  and  quadratures  are  feedback  methods.  The  feedback 
methods  are  increasingly  in  the  focus  of  interest  in  the  area  of  numerical 
solution  of  PDE  and  computational  mechanics.  Many  reports,  preprints  and 
papers  dealing  with  the  adaptivity  appeared  recently.  It  seems  to  be  that 
successful  computation  of  today's  large  complex  problems  cannot  be  made  with¬ 
out  some  kind  of  feedback  or  adaptive  approach. 

In  general  we  call  a  numerical  process  a  feedback  process  when  it  uses 
the  data  in  a  sequential  way  so  that  the  flow  of  computation  is  determined  by 
some  feedback  from  the  previous  information  obtained  during  the  computation 
(see  eg.  [18]  [19]  [20]).  However,  such  a  feedback  procedure  is  advantageous 
only  if  it  has  some  useful  properties ,  for  example ,  if  it  produces  the  re¬ 
quired  information  by  less  work  and  higher  reliability  than  the  nonadaptive 
process.  In  order  to  study  the  adaptive  procedures,  more  precise  definition 
of  the  used  notions  has  to  be  introduced,  see  eg.  [1]  [3]  [16], 

It  is  worthwhile  to  distinguish  between  the  feedback  and  adaptive 
procedures.  A  feedback  procedure  (method)  is  said  to  be  adaptive  if  it  is 
optimal  with  respect  to  certain  performance  measure;  hence  the  adaptivity  is  a 
relative  notion. 

The  analysis  of  a  method,  adaptive  or  nonadaptive,  needs  an  exact 
description  of  the  class  of  problems  which  the  method  is  designed  for.  The 
comparison  of  a  feedback  and  nonfeedback  approach  is  also  relative  to  the  set 
of  the  problems  under  consideration.  It  has  been  shown  that  if  the  class  of 
the  problems  is  big  enough,  then  in  the  worst  case  for  the  class  (or  even  in 
the  average  case)  the  non-feedback  methods  are  as  good  as  the  adaptive 
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methods.  See  e.g.  [18]  [20].  In  this  case  one  can  conclude  that  the  non¬ 
feedback  method  will  be  computationally  more  efficient  than  the  feedback  one 
(especially  for  parallel  computation).  On  the  other  hand,  as  experience 
shows,  for  the  majority  of  practically  important  classes  of  problems  adaptive 
methods  perform  much  better,  they  are  more  effective  and  reliable  than  the 
nonadaptive  ones. 

2.  FEEDBACK  AND  ADAPTIVITY 

We  discuss  now  the  notion  of  feedback  and  adaptivity  in  numerical  compu¬ 
tation  in  more  detail.  These  concepts  were  discussed  in  connection  with  the 
finite  element  method  in  [1]  [9].  For  a  more  general  concept  see  [16]. 

We  will  illustrate  the  introduced  notions  on  simple  example  of  the  two 
point  boundary  value  problem: 

(2.1a)  -(au')"  +  bu  *  f 

(2.1b)  u(0)  -  u(l)  -  0 

where  the  aim  is  to  find  the  approximate  solution  with  a  sufficient  accuracy 
measured  in  the  energy  norm. 

We  will  illustrate  the  introduced  notions  in  the  connection  with  the  h- 
version  of  the  finite  element  method  (with  elements  degree  p  »  1). 

In  general,  let  P  be  the  class  of  problems  to  be  solved.  (In  our  case 
this  class  is  characterized  by  class  of  functions  (a,  b,  f) .) 

Denote  by  S  the  set  of  numerical  solutions .  (The  set  of  all  piecewise 
linear  functions  of  Hq.) 

An  algorithm  is  a  mapping  which  maps  the  set  of  problems  P  into  the  set 
S  .  More  precisely,  an  algorithm  relates  to  an  information  operator  N  which 
may  be  referred  to  a  set  of  subroutines  such  that  the  information  needed  by 


the  algorithm  may  be  computed  for  the  given  problem  P  (e.g.  the  evaluation 
of  a,  b,  f) .  It  should  be  noticed  that  the  information  is,  in  general,  only 
partial*  (Values  of  a,  b,  f  can  be  computed  at  certain  points  and  with 
certain  precision  only.) 

There  are  two  kinds  of  information  operator  (see  [18]  [20]);  the  first  is 
the  non-feedback  one  (in  [18]  it  is  called  non-adaptive  information),  which 
uses  no  information  of  computed  results  of  the  solution  (e.g.  when  using  a 
uniform  mesh).  The  other  is  the  feedback  information  operator  which  uses 
sequentially  the  information  obtained  earlier  (e.g.  the  solution  on  the  pre¬ 
viously  used  meshes).  However,  no  matter  what  kind  of  Information  is  used,  an 
algorithm  produces  a  sequence  of  states  (the  sequence  of  meshes)  and  a 
sequence  of  corresponding  solutions  (the  finite  element  solutions).  The  pairs 
of  the  states  and  solutions  compose  the  traj ectory  of  the  algorithm.  Often  we 
also  say  that  the  sequence  of  the  states  is  the  trajectory.  The  map  mapping 
the  previous  states  into  a  new  state  is  called  a  transition  operator  construc- 
tion  of  the  new  mesh  from  the  old  ones) .  In  the  non-feedback  process  the 
transition  operator  of  course  makes  no  use  of  the  previous  states.  The  algo¬ 
rithm  is  completely  described  by  the  structure  of  the  transition  operator. 

In  summary  a  numerical  procedure  consists  of  the  following  components: 

1)  A  set  of  problems  P  ((1)  with  the  class  of  admissible  a,  b,  f). 

2)  The  set  S  of  the  numerical  solutions  (piecewise  linear  functions  of 


3)  The  set  of  states  X  (the  set  of  all  permissible  meshes  A). 

4)  A  set  of  information  I  (the  pointwise  values  of  a,  b,  f)  and  an 
information  operator  N  related  to  the  algorithm  N:  H  *  P  x  X  *  S  +  I. 

(The  subroutines  which  obtain  the  information  of  the  problem  p  (  P  at  n-th 
stage,  n  €  N,  using  the  meshes  A^  €  X  and  FE  solution  u^  €  S»  i  <  n-1). 


4 


5)  Transition  operator  T,  T:  N  *  I  ♦  X.  (Using  information  I  the 
transition  operator  creates  a  new  mesh.) 

6)  Solution  operator  I,  R:  X  x  I  +  S.  (Obtaining  the  FE  solution  for 
the  new  mesh  using  information  of  the  problem.) 

7)  A  set  of  trajectories  T  which  consists  of  all  sequences 

w  -  ((Aq.Uq)*  (Aj tu  ),  ...  ). 

An  algorithm  maps  any  element  P  (  P  into  a  trajectory  u>  G  T,  We  mention 
that  we  did  not  include  (for  simplicity)  the  stopping  criterion  in  the 
process . 

A  feedback  algorithm  is  useful  only  if  it  performs  better  than  any  non¬ 
feedback  one  for  the  set  P.  In  [1]  and  [14]  the  notion  of  a  performance 
measure  was  introduced.  A  performance  measure  p  is  a  map  from  the  set  of 
trajectories  into  R1. 

A  feedback  algorithm  is  said  to  be  adaptive  (relatively  to  P, 

S,  X,  u)  if  it  is  optimal  for  the  performance  measure  p.  (In  general  such 
optimal  algorithm  doets  not  have  to  exist  for  some  particular  (P,  S,  X,  p). 

Several  concrete  performance  measures  were  introduced  in  [l]  for  the 
finite  element  method.  The  convergence  measure  p  was  defined  as  binary 
measure  p:  I  ♦  {0,1}.  If  the  solution  converges  then  p(a>)  *  1,  otherwise 
p(oj)  «  0.  A  convergence  rate  measure  was  defined  as  follows:  For  any 
problem  P  €  P  let  <&(n,P)  be  the  minimal  error  obtainable  by  the  mesh  A 
with  m  +  1  nodal  points,  i.e.,  dim(S)  *  m  where  S  is  the  finite  element 
space  for  given  mesh.  The  particular  feedback  method  is  optimal  p(u>)  “1  if 

e(A,,P) 

(2.2)  ♦(m(A.),P)  *  K(P)’  J  " 
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i.e.,  when  the  error  e  of  the  finite  element  solution  on  the  sequence  of 
meshes  A^  constructed  by  the  algorithm  is  of  the  same  magnitude  as  the  theo¬ 
retically  smallest  possible  error  4  for  all  P  €  P.  In  [9]  we  described 
some  concrete  algorithms  which  are  adaptive  with  respect  to  the  convergence 
measure  (defined  above)  for  the  set  Pj  consisting  of  all  problems  with 
0  <  Sq  <  a(x)  <  a  <  -,  0  <  b(x)  <  b^  <  ®  and  f  such  that  the  exact 

solution  belongs  to  H*. 

These  algorithms  are  adaptive  with  respect  to  the  convergence  rate  only 
if  the  set  P^  is  restricted  to  a  smaller  set  This  set  ?2  never¬ 

theless  broad  enough  to  Include  vast  majority  of  problems  with  practical  sig¬ 
nificance.  In  [9]  we  have  shown  an  example  of  P  €  Pj,  p  ^  ^2  suc^  that  the 
algorithm  produces  meshes  with  (2.2)  unbounded.  Algorithm  of  a  similar  type 
as  discussed  in  [9]  was  Implemented  by  the  code  FEARS  [15]  and  PLTMG  [10]  in 
the  two  dimensional  setting. 

In  [14]  we  theoretically  analyzed  the  feedback  and  adaptivity  of  some 
algorithm  based  on  the  h-p  version.  We  have  proven  the  adaptivity  with  re¬ 
spect  to  the  convergence  measure  for  a  large  class  P.  The  convergence  rate 
measure  was  defined  slightly  different  than  before.  We  considered  the  set 
Pj  of  problems  with  solutions  which  are  analytic  except  a  finite  number  of 
points  Xq  where  the  solution  is  of  the  type  |x-  Xq|°,  o>  .5. 

We  have  shown  that  for  any  P  f  P^  and  any  mesh  A  and 
p-distribution 

(2.3)  4(n,P)  >  Const(P)  r^0”*^11 

n 

where 

rQ  -  (/T  -  l)2 

and  n  is  the  number  of  degrees  of  freedom  (i.e.,  dimension  of  the  used 
finite  element  space). 


The  algorithm  which  we  studied  and  implemented  in  [14]  leads  to  an  error 
2.4)  e  *  const(P) ($(n,P))Y 

rith  y  independent  of  P  and  a.  Combining  (2.3)  and  (2.4)  we  immediately 
lee  that  our  algorithm  is  adaptive  with  respect  to  and  obviously  defined 

:he  convergence  rate  measure.  We  can  also  readily  see  that  any  non-feedback 
ilgorithm  cannot  give  an  exponential  rate  of  convergence  for  an  arbitrary 
jroblem  in  this  set. 

Various  feedback  algorithms  appeared  in  the  literature  and  are  available 
is  codes.  Nevertheless  there  is  no  available  analysis  of  the  adaptive  features 
in  the  sense  explained  above  although  numerical  experimentation  with  these 
codes  can  give  very  good  insight  and  lead  to  the  proper  conjectures. 

3.  A  PARTICULAR  MODEL  PROBLEM 

Let  us  be  interested  in  the  nonlinear  problem 

(3.1a)  ~a(u',u,x)  +  b(u,x)  «  f(x,A)  x  €  I  =*  (0.1) 

(3.1b)  u(0)  -  u(l)  -  0 

with  the  aim  to  find  the  solution  u(x,A)  for  x  €  I  in  the  given  accuracy 
t  in  the  H*  norm  for  all  0  <  A  <  A^.  We  will  describe  briefly  the  main 
parts  of  a  feedback  procedure  for  solving  (3.1a)  (3.1b).  Although  (3.1)  is 
one  dimensional  problem,  the  main  ideas  which  will  be  explained  here  are  not 
one  dimensional  and  they  are  being  implemented  in  a  two  dimensional  setting. 
The  solution  is  computed  by  the  continuation  method  using  (for  a  fixed  mesh) 
adaptive  procedures  as  PITCON  (see  eg.  [17]).  The  problem  (3.1)  has  some 
common  features  with  the  methods  for  solving  parabolic  equations.  The  con¬ 
tinuation  parameter  g  plays  analogous  role  as  the  time  t  in  the  parabolic 
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oblem  and  will  be  called  "time"  for  convenience.  Of  course  there  are  also 
sential  differences  here  too.  For  more  about  adaptive  procedures  for  para- 
ilic  equations  we  refer  e.g.  to  [11]. 

In  the  following  we  shall  discuss  a  feedback  mesh  regridding  procedure, 
is  clear  that  while  the  "time"  parameter  X  advances  a  fixed  mesh  in 
sneral  cannot  work  well  and  the  mesh  is  necessary  to  be  regrldded.  However, 

10  frequent  change  of  mesh  also  increases  the  cost  of  computation.  Therefore 
s  shall  study  the  regridding  strategy  (including  the  criterion  of  mesh 
tality)  as  well  as  the  adaptivity  in  this  procedure. 

THE  MESH  DENSITY  AND  OPTIMAL  MESH 

We  now  discuss  on  the  optimal  mesh  for  the  linearized  problems.  We 
isurae  that  the  h-version  FEM  Is  used  for  solving  the  problem,  and  the  poly- 
>mlal  degree  of  the  elements  is  assumed  to  be  p.  For  one  dimensional  case 
lis  problem  was  considered  in  [7],  [13],  by  using  the  graded  meshes. 

A  graded  mesh  A  is  obtained  by  a  given  grading  function  g(x)  and 
itensity  m(A)  (i.e.,  number  of  elements  in  A).  For  example,  on  [0,1], 
le  nodal  points  of  a  mesh  A  are 

>•1)  XI  “  g(“) »  i  *  0,1 ,2, ... ,m 

lere  g  is  a  strictly  increasing  continuous  function  with  g(0)  “0,  g(l)  ■ 
,  The  optimal  grading  function  is  directly  related  to  some  derivatives  of 
\e  function  u(x)  (see  e.g.  [7]  [8]).  We  have  shown  in  [13]  (see  also  [8] 

>r  example),  that  if  u(x)  *  xa,  a  >  .5,  the  function 

g(x)  =  xS 


8 


p  +  .5 
a  -  .5 
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he  optimal  grading  function.  It  is  easy  to  show  that  in  this  case 


)  ■  c|u<p+1)(x)|2p+1 

•e  g  *(x)  is  the  inverse  function  to  g(x).  The  similar  expression  holds 
[eneral  under  proper  assumptions  (see  [7]  {8]). 

This  approach  is  obviously  only  one  dimensional.  Let  us  modify  it  so 
:  it  can  be  easily  transfered  into  more  dimensions.  For  this  reason  we 
L  introduce  the  notion  of  the  mesh  density  function  (briefly  the  density) . 

Let  J  c  I  be  any  interval  and  {A^}  be  a  sequence  of  meshes.  Let 
^,J)  be  the  number  of  elements  of  the  mesh  Ak  contained  in  J  and 
jme  that 


3) 


6(J) 


lim 

k-*° 


m(Ak,J) 

"(V 


well  defined. 

Function  6(J)  can  be  extended  to  a  a-additive  set  function  and  we  will 
ume  that  it  is  differentiable.  Denote  by  S(x)  its  derivative.  This 
ns  that 


4) 


<5(x) 


lim 
x€  J 


j|+0 


6(J) 

W~ 


sts  for  any  x  €  I.  We  will  call  the  function  S(x)  the  mesh  density 
ction  for  the  sequence  of  meshes  A^.  If  the  grading  function  exists,  then 
is  easy  to  see  that 

5)  <$(x)  -  g  *(x) 


hence 


Clu<P+1,U) 


2 

|  2p+l 


6) 


6(x) 
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TABLE  1. 


A 

m 

e 

X 

mesh 

action 

P2(k) 

1 

3.0 

64 

0.04003 

1.3599 

1.1176 

not  O.K. 

regridded 

1 

3.0 

48 

0.04458 

1.0002 

1.0591 

O.K. 

3.7382 

2 

3.5 

48 

0.04888 

1.0293 

1.0114 

O.K 

2.4349 

3 

4.0 

48 

0.05638 

1.1069 

0.9417 

not  O.K 

regridded 

3 

4.0 

55 

0.04086 

1.0058 

1.1063 

O.K. 

3.1482 

4 

4.5 

55 

0.04222 

1.0120 

1.0882 

O.K. 

2.6894 

5 

5.0 

55 

0.04682 

1.0746 

1.0335 

O.K. 

2.4239 

6 

5.5 

55 

0.05469 

1.1803 

0.9562 

not  O.K. 

regridded 

6 

5.5 

58 

0.03977 

1.0180 

1.1213 

O.K. 

2.8321 

7 

6.0 

58 

0.03821 

1.0023 

1.1440 

O.K. 

2.6480 

8 

6.5 

58 

0.03955 

1.0386 

1.1244 

O.K. 

2.5254 

9 

7.0 

58 

0.04330 

1.1188 

1.0746 

not  O.K. 

intact 

2.4417 

10 

7.5 

58 

0.04900 

1.2320 

1.0101 

not  O.K. 

intact 

2.3844 

11 

8.0 

58 

0.05638 

1.3689 

0.9417 

not  O.K. 

regridded 

11 

8.0 

57 

0.03963 

1.0673 

1.1232 

O.K. 

2.6380 

12 

8.5 

57 

0.03452 

1.0259 

1.2035 

not  O.K. 

intact 

2.5877 

13 

9.0 

57 

0.03099 

1.0052 

1.2703 

not  O.K. 

regridded 

13 

9.0 

44 

0.04460 

1.0076 

1.0589 

O.K. 

2.6540 

14 

9.5 

44 

0.04247 

1.0003 

1.0850 

O.K. 

2.5827 

15 

10.0 

44 

0.04083 

1.0089 

1.1066 

O.K. 

2.5254 

16 

11.0 

44 

0.03800 

1.0571 

1.0147 

O.K. 

2.4974 

•  •  • 

•  •  • 

•  •  • 

•  •  • 

T 


It  is  possible  that  the  new  mesh  could  be  unsuccessful,  for  example,  it 
Id  cause  divergence  in  the  iteration  of  finding  the  solution,  or  it  could 
e  a  solution  with  low  accuracy.  In  these  cases  one  has  to  adjust  the 
h.  We  will  not  go  into  detail  of  the  adjusting  procedure. 

According  to  our  assumption,  when  considering  a  process  which  constructs 
imal  mesh  at  each  step,  the  ratio  P2(k)  -  1  +0  «  3.  Our  feedback 
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particular  let 


1.2)  Dn(x,X)  -  - T 

u  1  +  10[x  -  sin  X/10] 


i  select  the  following  parameters  steering  the  described  feedback  process: 


error  tolerance 

T 

■ 

.05 

maximal  mesh  quality  indicators 

V.ax 

- 

1.2 

*max 

1.1 

the  cost  control  parameters 

*1 

=* 

2 

*2 

- 

1.2 

the  cost  exponent 

U 

- 

2 

the  cost  regridding  factor 

9 

3 

2 

the  maximal  number  of  predicted  steps 

kjnax 

- 

10 

ie  mesh  is  constructed  in  a  binery  tree  structure  to  simulate  data  structure 
led  in  the  two  dimensional  setting.  This  means  that  any  node  of  the  mesh  can 
>  obtained  by  a  successive  bisection  of  elements.  The  base  mesh  which  is 
xed  in  the  regridding  process  has  level  3  and  hence  it  consits  of  8 
iterval  of  the  size  1/8. 

The  results  are  shown  in  Table  1.  We  see  that  the  regridding  is  trigered 
’  various  reasons 

a)  The  error  is  not  acceptable  (i|>  <  1):  see  mesh  No.  3,  6,  11. 

b)  The  mesh  shape  is  bad  (x  >  Xmax)!  see  mesh  No.  1. 

c)  The  mesh  intensity  is  too  high  (ij>  >  <J»max) :  see  mesh  No.  13. 

Sometimes  the  mesh  quality  indicators  exceed  the  given  criteria  but  the 
>st  analysis  shows  that  the  mesh  is  better  being  kept  intact,  see  mesh  No.  9, 
),  12. 


ntensity  of  the  optimal  mesh  at  X  (to  achieve  the  given  accuracy  t) 


et  Xj,  j*l,2,...  be  the  continuation  steps  produced  by  PITCON,  for 
xample.  The  ideal  total  cost  (in  k  steps)  will  be 


10.4) 


k 


l  C(o  (X  )) 
j-0  U,T  2 


We  will  consider  the  ratio 


10.5) 


P2(k) 


CCloCal(k)  ' 


ind  the  measure  1  if  P2(k)  <  Y  (for  some  given  X  >  1)  for  all 

■  >  kg,  p2  “  0  otherwise. 

We  can  also  consider  the  measure  with  respect  to  the  time  inverval 
^1*^2^  an  ana^°8ue  ^y. 

The  important  question  is  whether  P2(k)  <1+0,  i.e.,  the  cost  of  our 
feedback  is  less  than  the  cost  when  the  mesh  is  changed  to  be  optimal  at  each 
itep. 


The  adaptivity  of  the  proposed  feedback  can  be  analyzed  under  various 
3tringent  assumptions.  Nevertheless  we  will  not  go  into  this  analysis. 
Instead,  we  show  a  numerical  example  computed  by  a  code  we  wrote. 


LI.  A  NUMERICAL  EXAMPLE 

As  we  have  seen,  the  entire  process  is  characterized  by  the  (unnormal- 
Lzed)  density  function  Dq(x,X).  We  can  simulate  the  process  of  regridding 
assuming  that  the  degree  of  element  p  and  Dq(x,X)  is  given.  Then  for  p  » 
l  the  error  indicator  n(I)  is  given  by 


we  have  to  define  exactly  the 


To  define  properly  the  measure  Up 
cost.  There  are  various  possibilities.  We  will  adopt  the  following 
definitions : 


Let  Xj,  j  =  1,2,...  be  the  continuation  steps,  Aj  be  the  meshes  at 

Xj.  These  meshes  are  not  necessarily  different,  we  assume  that  the  regrid- 

dings  only  occur  at  X.  ,  v  -  1,2,...  .  Then  we  define  the  total  cost  in 

Jv 

k  continuation  steps  to  be 


(10.1) 


r 

total 


(k) 


k 

I  (1  +  a  0)C(m  ) 
j-1  3  3 


where  Oj  “  1  if 
and  nij  ■  m(Aj) 

We  can  also 


j  ■  Jv  (i.e.,  regridding  occurs),  ■  0  otherwise 
is  the  intensity  of  mesh  Aj . 

consider  the  cost  in  certain  time  interval  [A^.A^]: 


) 


(10.2) 


Ctotal(Al»A2* 


l  (1  +  a, e)C(m, ). 

Al<Xj<A2  3  3 


The  cost  of  the  feedback  and  non-feedback  methods  will  be  denoted  by 
superscripts  a  and  n  respectively,  then  we  define 


(10.3) 


Pj(k) 


Ctotal 

cn 

utotal 


(k) 

(k) 


and  “  1  if  P^(k)  <  y  (for  some  given  0  <  y  <  1)  for  all  k  >  kg, 

”  0  otherwise. 

An  analogue  of  (10.3)  may  also  be  made  for  the  total  cost  defined  by 

(10.2). 

Another  measure  can  compare  the  cost  of  the  introduced  feedback 

procedure  with  the  ideal  cost  of  the  method  when  optimal  meshes  are  used  at 
every  step  (neglecting  the  cost  of  regridding).  Let  m  (X)  be  the 

U  -  T 
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10.  THE  ADAPTIVITY  OF  THE  REGRIDDING  PROCESS 

It  is  clear  that  the  process  described  in  Section  9  is  a  feedback 
process.  We  will  first  include  it  into  the  frame  discussed  in  Section  2.  The 
set  of  problems  is  defined  by  the  set  of  functions  a,  b,  f  in  (3.1a).  The 
solution  set  is  the  finite  element  space  of  elements  of  degree  p  on  the 
meshes  A.  The  states  are  the  pairs  (X,A)  where  A  is  the  mesh  at  time 
X.  The  information  operator  consists  of  a  set  of  subroutines  which  evaluate 
the  Involved  functions  to  obtain  necessary  information  described  in  the  last 
section.  The  transition  operator  makes  decision  of  regridding  and  construct 
new  meshes  when  necessary.  The  solution  operator  then  computes  the  finite 
element  solution  on  the  given  mesh  X. 

In  solving  the  nonlinear  problem  by  continuation  (for  example,  by 
PITCON) ,  another  feedback  process  is  involved.  We  shall  not  address  here  the 
question  of  adaptivity  of  this  continuation  process. 

The  regridding  process  itself  is  very  complex.  For  the  adaptivity 
analysis  various  performance  measures  can  be  considered.  Numerical  exper- 
iments  can  be  useful  for  formulation  some  conjectures  for  theoretical 
investigations . 

Like  the  convergence  measure  discusssed  earlier,  we  introduce  here  a 
natural  performance  measure  Pq  such  that  if  all  solutions  in  the  trajectory 
satisfy  given  tolerance  r  (which  may  be  restricted  to  some  range,  say  1Z- 
10Z) ,  then  Pq  ■  1,  otherwise  Uq  ■  0.  The  algorithm  is  said  to  be  adaptive 
with  respect  to  pg  if  Pq  *  1  for  all  problems  P  €  Pq.  It  is  clear  that 
our  algorithm  is  adaptive  for  a  large  class  of  problems  Pq. 

Another  natural  performance  measure  compares  the  cost  of  our  feedback 
procedure  with  the  cost  of  the  computation  when  fixed  mesh  is  used.  We  define 
the  measure  p^  so  that  p^  ■  1  if  the  cost  of  the  feedback  is  not  larger 
than  the  cost  of  the  fixed  mesh  approach.  Otherwise  p^  *  0. 


-  it  • 
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(8.4) 


ci  “  c(°0,t) 


where  qq  ^  is  Che  intensity  of  optimal  mesh  at  Aq. 

Assuming  that  the  average  step  AX  for  X  (which  is  selected  in  a 
feedback  fashion)  will  stay  the  same  as  in  the  past,  we  will  compute  a 
(X0*Xl)(x) 

j - 1  •>  and  t^e  corresponding  intensity  m  as  in  the  last 


density  D 


section  for  each  A  «  A^  +  kAA,  k  ■  0,1,2,...  .  The  predicted  cost  for  this 
mesh  is  then 


,k  _  <VX>- 


(X  -  XQ  +  kAX). 


This  procedure  is  proceeded  until  k  ■  k' : 


Ck  >  min{C  /k,  k.C  }, 

P  C  L  1 


1  Ci  ’ 


and  >  1  are  given  a-priori,  C°  is  the  ideal  cost  when  last  mesh 

regrldding  occurred.  The  predicted  mesh  will  be  constructed  from  the  density 
(An,A)  (XX) 

D  (x)  and  m  with  A  ■  A.  +  (k'  -  1)AA, 

P  U 

Roughly  speaking,  the  meaning  of  the  above  selection  is  that  the  cost  of 
a  new  mesh  should  not  exceed  too  much  the  cost  of  the  ideal  mesh  and  should 
take  into  account  the  change  of  the  optimal  mesh. 


9.  THE  MAIN  FLOW  OF  REGRIDDING  PROCESS 

Here  we  briefly  show  the  main  flow  of  the  feedback  algorithm  for  the 
design  of  the  mesh. 
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6)  Make  a  cost  analysis  for  determining  A^.  (See  Section  8.) 

7)  Compute  the  uncertainty  function  (p(x,A)  of  the  prediction  by 
comparing  the  expected  density  with  the  ideal  one  (after  being  computed). 

8)  Design  the  new  mesh  if  the  current  mesh  has  to  be  changed. 


8.  THE  COST  ANALYSIS 

The  cost  of  the  adaptive  mesh  consists  of  two  parts: 

1)  The  cost  of  computation  by  the  continuation  method  (for  example  by 
PITCON)  per  step: 

(8.1)  C  -  C(m)  -  mU  . 


where  C(m)  is  the  cost  function  which  is  supposed  to  be  only  depending  on  the 
Intensity. 

2)  The  cost  of  regridding  the  mesh  which  will  be  assumed  as  a  multiple 
of  the  cost  of  the  solution: 

(8.2)  CQ  -  0C(m) . 


A  mesh  is  considered  to  be  changed  if  either  it  cannot  give  the  required 
accuracy  or  it  costs  too  much  (although  the  error  is  admissible).  The  first 
case  is  characterized  by  <  1,  and  the  latter  is  determined  by  the 
criteria  x  >  Xmax  or  x  >  <J»oax  when  x^  and  *max  are  given  a-priori. 

We  define  the  cost  for  the  current  mesh  (at  Aq) 

I®  if  X  <  1, 


(8.3) 


C 

c 


|C(mc)  otherwise. 


where  mc  is  the  intensity  of  the  current  mesh. 

Furthermore,  the  ideal  cost  (at  Aq)  is  given  by 


.  *  .  *  .  *  '  .  ’ 
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turbance  of  an  optimal  mesh  causes  changes  of  accuracy  of  second  degree  (see 
[7]),  we  do  not  need  to  know  the  optimal  density  too  accurately.  It  is 
computationally  advantageous  to  work  only  with  meshes  which  have  density  given 
on  the  base  mesh.  For  example,  if  the  elements  of  degree  p  ■  2  are  used  for 
the  computation  of  the  solution,  we  will  assume  that  the  density  is  given  by 
the  same  2  degree  element  function  on  the  base  mesh  and  the  density  is  given 
by  its  values  in  the  nodal  points  of  the  base  mesh. 


The  prediction  of  the  density  goes  as  follows : 

1)  Define  piecewise  constant  function  Dq(x,A), 


X  <  Aq  on  the  mesh 


x  f  J 


2)  By  a  fitting  procedure  construct  the  (unnorraalized)  density  function 
Dq(x,A),  called  ideal  density,  on  the  base  mesh. 

3)  Given  function  Dq(x,A)  for  A  <  A_  (which  is  given  by  its  values 
in  the  nodal  points  of  the  base  mesh) ,  extrapolate  it  for  A  >  A^  (by  an 
extrapolation  technique  which  preserves,  for  example,  its  positivity),  and 
denote  it  by  Dg(x,X)  (called  the  expected  density). 

4)  A  new  mesh,  called  predicted  mesh,  is  designed  for  A^  <  A  <  A^,  by 

(6.2) 

(A0,X  ) 

D  (x)  *  max  D  (x,A). 

P  xQ<x<x1  e 

(VV 

5)  Compute  the  Intensity  m  leading  to  the  prescribed  accuracy 

(Xq  ,Xj) 

using  mesh  D  (x) 
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and 


(6.3) 


«0.Ji> 


(  x) 


(X  ,x  ) 

°o  <*> 

<VV  ' 


Ur 


Furthermore,  the  intensity  used  in  [Ag,A^I  will  be  determined  by 


u0»xi) 

m 


max 


l0+W2prr  (5o(«.x))2^  >«p 


(/ 


Vuxi  t1/p  1  (s 


a°'Al,(x»2P 


6*) 


7.  DENSITY  PREDICTION 

The  problem  (3.1)  is  solved  by  a  continuation  method  with  respect  to  the 
parameter  X.  This  means  that  we  know  the  approximate  solution  for  A  <  A^ 
and  want  to  predict  the  behaviour  of  the  solution  for  A  >  A^  and  design  an 
optimal  mesh  for  XQ  <  X  <  X^,  where  Aj  will  be  properly  chosen.  The 
prediction  process  is  essentially  a  pattern  recognition  which  is  composed  by 
learning,  classification,  prediction  and  correction.  Some  ideas  of  pattern 
recognition  were  used  in  [11]  where  the  concepts  of  the  shape  and  intensity  of 
a  mesh  were  introduced. 

As  we  have  seen  above  the  shape  of  the  mesh  is  determined  by  its 
density.  The  function  Dq(x,A)  defined  by  (A. 10)  is  an  optimal  density  up  to 
a  multiplicative  factor  — .  It  also  contains  the  information  about  optimal 
intensity  using  (5.5)  with  its  approximate  version 


(J 


1 

3p 


We  will  approximate  Dg(x)  by  a  (finite  element)  function  on  a  relatively 
coarse  mesh  called  base  mesh  which  does  not  change  with  X.  Because  a  per- 
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Using  the  formulae  of  Section  4  we  get 


(5.3) 


2p  2 

m  rT 
T 


2p+l  { 

“o  I 


(y*>) 


2p+l 


I  6(x) 


2p 


dx 


(5.4) 


2p  2 

m  re 
c 


*^i  ,  (v*»2p+1 


/ 

I  <S(x) 


2p 


dx 


(5.5) 

Hence  we  have 

(5.6) 

and 

(5.7) 


2p  2  2p+l 

0,t  0 


*  •  (|) 


[■ 


1/P 


(  I  h  <j» 

2p+l 


12P+1  1/, 


1/2 


x  .  a  (eQ(,));  -  axfp 

I  (6(x))  2p 


6.  OPTIMAL  MESH  FDR  SIMULTANEOUS  APPROXIMATION 

In  previous  sections  we  discussed  the  problem  of  the  optimal  mesh  for 
approximating  u(x,A)  for  one  fixed  value  of  A.  We  would  like  to  get  an 
optimal  mesh  for  set  of  functions  u(x,A),  A  *  (Aq.A^,  To  this  end  we  could 
use  (4.14)  and  minimize 


(6.1) 

among  all  densities 
mately  we  will  take 


max 

Aq<A<A 

<5(x)  . 


I 

L  1 

To  find  such 


(s  0u,i»2'*1 

— - - - -  dx 

(6(x))2p 

6(x)  is  not  easy. 


Hence  approxi- 


(Xq,A  ) 

Dq  1  (x)  -  max  Dq(x,A) 


(6.2) 
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(4.15) 


500^ 


*0  *  ^0 
1  +  <f>o 


For  large  uncertainty  <p  q  »  1  we  get  by  this  technique  a  nearly  uniform 
mesh. 

Finally  let  us  mention  that  under  sufficiently  strict  assumptions  our 
heuristic  reasoning  could  be  given  a  rigorous  base.  In  addition,  we  underline 
that  although  we  explained  the  approach  in  one  dimension,  it  is  valid  in  more 
dimensions,  too. 


5.  THE  MESH  QUALITY  INDICATORS 

In  the  feedback  process  we  would  like  to  use  one  mesh  for  largest  pos¬ 
sible  range  of  \  because  of  the  cost  consideration.  This  leads  to  the  need 
of  quantitative  expression  of  the  quality  of  the  mesh.  We  will  characterize 
the  quality  by  two  indicators.  The  first  one  is  related  to  the  intensity  and 
second  to  the  density. 

First,  we  let 

(5.0  *  • 

T 

where  mc  is  the  intensity  of  the  current  mesh  and  mT  the  intensity  of  the 
mesh,  having  the  same  density,  achieving  (exactly)  e(A)  ■  x.  For  determining 
m  (4.14)  is  used.  If  e(A)  >  t,  then  the  mesh  is  not  acceptable  and  we 

T 

have  ip  <  1.  If  ip  »  1,  then  the  used  mesh  is  too  fine  and  the  cost  can  be 
too  high. 

We  define  the  second  indicator  by 


where  Oq  T  is  the  intensity  for  the  optimal  mesh  and  accuracy  t.  Clearly 
we  have  x  *  1»  and  x  »  1  indicates  that  the  current  mesh  has  a  density 
which  is  too  far  to  be  optimal. 
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the  elememts  of  the  mesh  A^.  Then 


(4.11) 


60(x) 


60(x) 


Is  the  density  of  the  optimal  mesh  where 

(4.12) 


-0  "  /  Do(x)  dx* 


Dq(x)  will  be  called  the  unnormallzed  optimal  density. 

2 

Let  us  assume  that  the  global  error  e  (A)  can  be  written  in  the  form 

(4.13)  e2(A)  -  l  n2(J). 

J^A 

Then  for  the  density  6(x)  and  the  intensity  m  *  m(A)  we  have  the  relation 

I'.  ,  . >2p+l 

,,  ,,x  .  .  2o  2,.,  2 p+1  <•  ^VX;j 

(4.14)  lim  m  *e  (A)  ■  uft  J  - z -  dx. 

m+<»  0  I  5(x)2p 

(4.14)  expresses  the  error  for  the  mesh  A  with  the  density  5(x)  through 
the  optimal  density.  Obviously  the  choice  S(x)  «  6q(x)  minimizes  the 
expression  on  the  right  hand  side  of  (4.14)  and  gives  the  minimal  value  ui2p+^ 
(because  /  6Q(x)dx  -1).  In  practice  we  express  6qCx)  by  the  error  indica¬ 
tors  and  (4.10)  (4.11).  Nevertheless,  the  error  indicators  (especially  in 
more  dimensions)  are  showing  a  dispersion  and  a  smoothing  technique  has  to  be 
used.  In  addition  it  is  useful  to  assume  that  we  have  in  (4.14)  for  disposi¬ 
tion  6q(x)  +  <p(x)  instead  6q(x)  °nly  where  <P(x)  expresses  the  uncer¬ 
tainty  in  the  determination  of  6q(x)  which  can  be  determined  by  the  smooth¬ 
ing  process  mentioned  above;  for  example,  we  can  assume  that  |<p(x)|  <  <Pq. 

Now  the  optimal  mesh  under  this  uncertainty  minimizes  (4.14)  for  the  worst 
case  of  the  uncertainty.  This  leads  to  the  density  600(x) 
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with  the  constant  C  such  that 


(4.7) 


/  6(x)  dx  -  1. 

I 


The  notion  of  the  density  function  can  be  obviously  easily  extended  to 
n-d  linens  ions . 

If  the  density  function  is  given,  then  it  is  not  hard  to  construct  a  mesh 
with  this  density  6(x)  and  given  intensity  m(A).  We  shall  assume  that  all 
meshes  we  are  dealing  with  have  a  density  function. 

Let  us  assume  now  that  an  error  indicator  n(J)  of  every  element 
J  €  A  of  the  mesh  is  computable  and  it  is  such  that 


(4.9) 


n(J)  -  >e"E(j)(1  +  o(1)) 


where  lellE(j)  ls  the  (local)  error  (in  the  energy  norm)  of  the  interpolant 
of  u  on  J.  Then  we  get 


(4.9) 


u*  -  C|u'^>U)| 

| j|+o  IjI2^1 


with  a-priori  known  constant  C.  See  [8 ]  for  more  details.  Similar  situation 
occurs  in  more  dimensions.  For  the  error  indicators  in  two  dimension  we  refer 
to  [2]. 

The  error  indicator  gives  through  (4.9)  the  information  about  the  deriva¬ 
tive  of  the  approximated  function  u  and  hence  also  the  density  for  the 
optimal  mesh.  To  this  end  we  define 


(4.10) 


Dn(x)  ■  lim 

0  |J|*0  I7 


1 _ 

Zp+1 


and  assume  that  Dg(x)  can  be  well  approximated  by  the  error  indicators  of 


regridding  process  obtains  the  ratio  p2<k)  *  2.5  (see  Table  1).  Thus  it  is 
more  efficient  that  the  previous  one. 

The  process  we  have  shown  is  of  course  not  the  only  possible  approach  and 
it  needs  detailed  theoretical  and  numerical  studies  in  the  spirit  we  explained 
above.  We  have  shown  it  as  an  example  of  application  of  principles  of  studying 
the  adaptivity. 

12.  CONCLUSIONS 

We  have  shown  the  main  questions  associated  to  the  design  and  analysis  of 
the  feedback  procedures.  The  theoretical  analysis  is  often  not  easy.  Never¬ 
theless,  every  feedback  algorithm  could  be  conjectured  to  be  adaptive  with 
respect  to  some  measures  and  class  of  problems  and  studied  theoretically  and 
experimentally  in  the  direction  we  explores. 
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The  Laboratory  for  Numerical  analysis  is  an  integral  part  of  the 
Institute  for  Physical  Science  and  Technology  of  the  University  of  Maryland, 
under  the  general  administration  of  the  Director,  Institute  for  Physical 
Science  and  Technology.  It  has  the  following  goals: 

•  To  conduct  research  in  the  mathematical  theory  and  computational 
implementation  of  numerical  analysis  and  related  topics,  with  emphasis 
on  the  numerical  treatment  of  linear  and  nonlinear  differential  equa¬ 
tions  and  problems  in  linear  and  nonlinear  algebra. 

•  To  help  bridge  gaps  between  computational  directions  in  engineering, 
physics,  etc.,  and  those  in  the  mathematical  community. 

•  To  provide  a  limited  consulting  service  in  all  areas  of  numerical 
mathematics  to  the  University  as  a  whole,  and  also  to  government 
agencies  and  industries  in  the  State  of  Maryland  and  the  Washington 
Metropolitan  area. 

•  To  assist  with  the  education  of  numerical  analysts,  especially  at  the 
postdoctoral  level,  in  conjunction  with  the  Interdisciplinary  Applied 
Mathematics  Program  and  the  programs  of  the  Mathematics  and  Computer 
Science  Departments.  This  includes  active  collaboration  with  govern¬ 
ment  agencies  such  as  the  National  Bureau  of  Standards. 

•  To  be  an  international  center  of  study  and  research  for  foreign 
students  in  numerical  mathematics  who  are  supported  by  foreign  govern¬ 
ments  or  exchange  agencies  (Fulbright,  etc.) 

Further  information  may  be  obtained  from  Professor  I.  BabuSka,  Chairman, 
Laboratory  for  Numerical  Analysis,  Institute  for  Physical  Science  and 
Technology,  University  of  Maryland,  College  Park,  Maryland  20742. 


END 

FILMED 

9-85 

DTIC 


